home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
The World of Computer Software
/
The World of Computer Software.iso
/
yamp16.zip
/
TEST.CPP
< prev
next >
Wrap
C/C++ Source or Header
|
1993-01-11
|
6KB
|
225 lines
/*
YAMP - Yet Another Matrix Program
Version: 1.6
Author: Mark Von Tress, Ph.D.
Date: 01/11/93
Copyright(c) Mark Von Tress 1993
Portions of this code are (c) 1991 by Allen I. Holub and are used by
permission.
DISCLAIMER: THIS PROGRAM IS PROVIDED AS IS, WITHOUT ANY
WARRANTY, EXPRESSED OR IMPLIED, INCLUDING BUT NOT LIMITED
TO FITNESS FOR A PARTICULAR PURPOSE. THE AUTHOR DISCLAIMS
ALL LIABILITY FOR DIRECT OR CONSEQUENTIAL DAMAGES RESULTING
FROM USE OF THIS PROGRAM.
*/
// make sure to define the global define variable IN_RAM if you want
// to use the medium model. This can be done in the compiler code
// generation window of compiler options. It may also be defined for
// the command line compiler
#include "virt.h"
/* required global declaration for the matrix stack object */
//unsigned int _stklen = STACKLENGTH;
MStack *Dispatch = new MStack;
VMatrix& function1(VMatrix &a, VMatrix &b)
{
// This function tests the freind functions and returns a value
a.Garbage("function1"); // check a and b
b.Garbage("function1");
Dispatch->Inclevel(); // increment push-pop level
VMatrix c("c", 1, 1); // create a local matrix
c = a + b + a + b; // check a repeated matrix addition
c.DisplayMat(); // print c
c = 431.2 + Tran(a) *b + 2.134; // check commutivity of scalar addition
c.DisplayMat(); // print c
c = -c; // check uniary minus
c.DisplayMat(); // print c
c = a - b; // matrix subraction
c.DisplayMat(); // print c
c = 5 - a - 5; // check commutivity of scalar subtraction
c.DisplayMat(); // print c
c = 5*a*5; // check commutivity of scalar multiplication
c.DisplayMat(); // print c
c = a % a; // check elementwise multiplication
c.DisplayMat(); // print c
c = a / 1234; // check scalar division
c.DisplayMat(); // print c
c = a / b; // check elementwise division
c.DisplayMat(); // print c
Dispatch->PrintStack(); // show stack before push
Dispatch->Push(c); // push c onto stack
Dispatch->PrintStack(); // examine stack after push
return Dispatch->DecReturn(); // decrement push-pop level, and return
// stack top
}
VMatrix &function0(void)
{
// test some of the output functions and raw matrix functions
Dispatch->Inclevel();
VMatrix d = VMatrix("d", 1, 1);
VMatrix a, c, H, hilb;
a = Reada("catchv.dat"); // read an ascii matrix
a.DisplayMat(); // display a
Writea("junk.dat", a); // write ascii matrix
a.Writeb("junk.bin", a); // write binary matrix
a = Readb("junk.bin"); // read binary matrix
a.InfoMat(); // display matrix info
a.DisplayMat(); // display a
d = Submat(a, a.r, 4, 1, 2); // take a submatrix of a
d.DisplayMat(); // display it
c = Ch(d, d); // horizontal concatenation
c.DisplayMat();
c = Cv(d, d); // vertical contatenation
c.DisplayMat();
c = Kron(Ident(3), d); // Kroniker's product
Setdec(1); // set number of decimals to print
Setwid(5); // set print width
c.DisplayMat();
H = Helm(4); // make a helmert matrix
H.DisplayMat(); // show that it is orthonormal
d = Tran(H) *H;
d.DisplayMat();
d = H*Tran(H);
d.DisplayMat();
a = Ident(4) + Fill(4, 4, 0.5); // redefine a
a.DisplayMat(); // print a
c = Tran(Helm(4)) *a*H; // diagonalize a
c.DisplayMat(); // display c
c = Inv(a); // invert a
c.DisplayMat(); // display c
VMatrix b; // create b
b = c*a; // redefine b
b.DisplayMat(); // display b
c = function1(a, a); // call a function
c.DisplayMat(); // display c
Dispatch->Push(c);
return Dispatch->DecReturn();
}
VMatrix ®ression(void) // do a multiple linear regression
{
Dispatch->Inclevel();
VMatrix a, xy, reg;
a = Reada("catchv.dat"); // read data
VMatrix b = VMatrix("b", a.r, a.c); //simplify indexes
int N = a.r; // N
int p = 3; // params
xy = Ch(Fill(N, 1, 1), Submat(a, N, 4, 1, 2)); // make x y
reg = Sweep(1, p, Tran(xy) *xy); // solve regression
// using sweep
reg.M(p + 1, p + 1) = reg.m(p + 1, p + 1) / ((double) (N - p));
// divide to get mse
VMatrix x = Submat(xy, N, 3, 1, 1);
VMatrix y = Submat(xy, N, 4, 1, 4);
VMatrix betahat = Inv(Tran(x) *x) *Tran(x) *y;
//solve regression using
//normal equations
betahat.DisplayMat();
reg.DisplayMat();
Dispatch->Push(reg); // put return mat on stack
return Dispatch->DecReturn(); // dec level & return
}
void altermatrix(VMatrix &t) // alter a matrix element
{
t.M(1, 1) = 2525.0;
VMatrix temp = MSort(t, 1);
temp.DisplayMat();
}
void testhuge(void)
{
//VMatrix VeryBig("vb", 1000, 1000);
cout << "Inverting " << endl;
VMatrix VeryBig = Inv( Helm(1000 ) );
VeryBig.InfoMat();
}
void version1p1(void)
{
//
// testreg.cpp for test of regression functions.
//
Dispatch->Inclevel();
VMatrix A = Fill(6, 6, 1.0);
VMatrix B = Fill(6, 1, 2.0);
(Vec(A)).DisplayMat();
(Vecdiag(A)).DisplayMat();
(Diag(B)).DisplayMat();
(Shape(A, 3)).DisplayMat();
(Sum(A, ROWS)).DisplayMat();
(Sumsq(A, COLUMNS)).DisplayMat();
(Cusum(A)).DisplayMat();
(Mmin(B)).DisplayMat();
(Mmax(B, ROWS)).DisplayMat();
VMatrix C = Ch(A, Diag(B));
Setwid(5);
Setdec(2);
C.DisplayMat();
Crow(C, 1, 0.6);
C.DisplayMat();
Srow(C, 1, 6);
C.DisplayMat();
Lrow(C, 2, 3,.5);
C.DisplayMat();
Ccol(C, 1, 0.6);
C.DisplayMat();
Scol(C, 1, 6);
C.DisplayMat();
Lcol(C, 2, 3,.5);
C.DisplayMat();
Dispatch->Cleanstack();
Dispatch->Declevel();
}
main()
{
VMatrix testit;
printf("%d\n", sizeof (testit));
testit = regression();
testit.DisplayMat();
VMatrix b = testit;
VMatrix a = 2*testit;
VMatrix c = a + b + a + b;
c.DisplayMat();
VMatrix t = function0() + function0();
t.DisplayMat();
altermatrix(testit);
testit.DisplayMat();
version1p1();
#ifndef IN_RAM
testhuge();
#endif
vclose();
}